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A NEW TWIST FOR THE SIMULATION OF HYBRID SYSTEMS USING THE 

TRUE JUMP METHOD 

ROMAIN VELTZ 


Abstract. The use of stochastic models, in effect piecewise deterministic Markov processes (PDMP), has 
become increasingly popular especially for the modeling of chemical reactions and cell biophysics. Yet, exact 
simulation methods, for the simulation of these models in evolving environments, ai‘e limited by the need to find the 
next jumping time at each recursion of the algorithm. Here, we report on a new general method to find this jumping 
time for the True Jump Method. It is based on an expression in terms of ordinary differential equations for which 
efficient numerical methods ai‘e available. As such, our new result makes it possible to study numerically stochastic 
models for which analytical formulas are not available thereby providing a way to approximate the state distribution 
for example. We conclude that the wide use of event detection schemes for the simulation of PDMPs should be 
strongly reconsidered. The only relevant remaining question being the efficiency of our method compared to the 
Fictitious Jump Method, question which is strongly case dependent. 

Key words. PDMP, piecewise deterministic Markov process; ODE, ordinary differential equations. 


1. Introduction. There is a growing number of processes (from biophysics, ecology, 
finance, internet traffic, queuing...) that are inaccurately modeled with ODEs. A classical 
example is the one of chemical reactions for which the reactants are well stirred but in such 
a small number that a continuous description is inappropriate. A description in term of a 
Markov process is more relevant in this case 0Gil77l . 

The study of Markov models, with intrinsic noise, comes with the computation of the 
distribution of the states at a given time. Other methods, aimed at understanding the noisy 
dynamics, rely on sample paths IIArn98l . Simulation methods which provide exact sample 
paths are very helpful for the computation of the state distribution. There are two broad 
classes of simulation methods 0GT13IICT97I . the Fictitious Jump (FJM) and the 

True Jump Method (TJM), see appendices. The last method was long known before being 
popularized by Gillespie 0Gil77l IGBOOl such that it bears the name of the Doob-Gillespie 
algorithm. 

There are several limitations to the above algorithms. The first occurs when the number 
of chemical reactions is large, the second occurs when there is a subpopulation with very fast 
kinetics compared to the other reactions and the third being if the transition rate^ are time 
dependent (in the case of the TJM). The first two limitations lead to very long simulation 
times while the last one requires to accurately solve an integral equation at each iteration 
which can be difficult and time consuming IIAnd07ll (but see IICT97I for a FJM version). 

The first two limitations can be addressed by using an approximate algorithm, for ex¬ 
ample the tau-leaping method IGilOll . or by using a time-scale separation IIKK13II which 
provides a PDMP as limit of the suitably scaled Markov model. In short, the limit is com¬ 
posed of an ODE, averaging the fast component of the Markov model and a jump process 
representing the slow component of the original Markov model. 

Finally, PDMPs IIDav84l IDav93l also arise in models in which one wants to couple de¬ 
terministic dynamics to a chemical process, as for example in the modeling of the mem¬ 
brane potential of a neuron coupled with the stochastic opening/closing of the ions channels 
IIBrel4l[PTW10l . The PDMPs are composed of two components: a continuous one described 
by the flow of an ODE and a jump process - that affects the continuous flow and describes the 
stochastic component of the process. Note that this class of Markov processes encompasses 
the one of chemical reactions with time dependent reaction rates as a subtype. 

* also known as rejection method 

^also called propensity functions in the case of chemical reactions 
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This motivates the need for the computation of exact sample paths of PDMPs. There are 
here two possibilities: 

1. If the continuous flow is known analytically, then the FJM is an exact method. Un¬ 
less the integral equation in the TJM can be solved analytically as well, the TJM is 
not exact in this case and should be discarded in favor of the FJM., 

2. If the continuous flow is not known analytically, none of both methods is exact. 
The error made in the simulation comes from the need to compute numerically the 
continuous flow, and in the case of the TJM, from the numerical computation of the 
jumping times. Hence, compared to the FJM, the TJM adds another error source. 
We shall see that our method adresses this last point. 

The simulation of PDMPs llACT+05l IBdSDlOl IGPB081 iRielJl is reviewed in IIRiel3l 
where the author justifies theoretically the original method of llACT~*~05l . albeit focusing only 
on the TJM. They propose a method based on an ODE solver with event detection IISha03l 
where the event is the integral equation that needs to be solved to And the jump instants (see 
appendix |^. Such numerical methods may not be available and suffers from two majors 
drawbacks. First, one may need to modify a library providing an ODE solver to implement 
the event detection. Given the complexity of such solvers IIHW96I . this may not be trivial. 
The second issue with the use of event detection scheme is that the time-step is only adapted 
to solve the ODE and not to locate precisely the event. This is a general word of caution 
concerning event detection for ODE, it is an ill-posed problem IISha03l and may lead to 
excessively large numerical errors as we shall see. Einally, note that all examples in IIRiel3l 
are non-stiff, hence being relatively easy to solve numerically. We address all these issues 
with one theoretical result. 

Here, we report on a new simulation algorithm derived from the TJM. It is based on 
the assumption that a regular ODE solver is available without further capabilities like event 
detection. The precision of our method is exactly the one of the ODE solver IISke86l IHW96I 
making the simulation of PDMP nothing more than the simulation of an ODE. In particular, 
global error estimation IISke86IIVis01IICP04l for ODE directly leads to global error estimation 
of the simulation of the PDMPs. Our new method is based on a change of time variable and 
recursively solving an ODE of dimension n-\- 1 where n is the dimension of the continuous 
flow of the PDMP. Hence, it comes with little additional cost given that one needs to simulate 
the continuous component of the PDMP anyway. 

The only relevant remaining question is how our new method compares to the EJM when 
the continuous flow is not known analytically. We will come back to this question and provide 
a numerical example. 

The paper is organized as follows. We first derive the theoretical result which is the basis 
of our method. Then we provide three numerical examples. In the first two examples, we 
show that the simulation methods based on event detection are to be reconsidered. We then 
provide a third example showing the effectiveness of our method compared to the EJM before 
drawing some conclusions. 


2. Description of the PDMPs. In this paper, we consider a PDMP by following IIGT13I 
in order to simplify the description of the process. Note that a more general definition is 
provided in IIDav84l|bav93l but we stick to flGT13l for simplicity. A PDMP is a jump process 
where the stochastic jump instants T„ are non decreasing, isolated with limit 
lim Tn = oo and the dynamics between the jumps is given by a flow of homeomorphisms 

n—>-oo 
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{(t>*)t>o on as follows: 

r Vf > 0, x{t) = E 

I ">0 

S 0 = To < Ti < Ta < • • • < oo, 

[ T„ < oo ^ T„ > T„_i and X(T„) ^ X(T„-). 

We assume that the continuous flow is produced by an ODE: 

( = F{M^)), t>0 

I Mx) = X. 


( 2 . 1 ) 


( 2 . 2 ) 


such that X F{x) is locally Lipschitz. The jump X{Tn) and the inter jump interval 
Tn — Tn-i are independent. The law of X{Tn) given X(Tn-i) = a; is n(a:, dy) and the law 
of the inter jump interval follows from 

P(T„ - T„_i > f|2f(T„_i) = x) = exp j Rtot{4>''[x))ds 

where the total transition rate function Rtot, is measurable and satisfies Rtot{x) > 0. We 
allow 


sup Rtot{x) < oo 

but Rtot must be locally bounded. This implies IIGT13I that there is a finite number of jumps 
on any bounded time interval. 

Remark 1. Using the usual trick, we can deal with the case where the vector field F in 
I is time-dependent, i.e. F : {x,i) —>■ F{x,t), but the existence of PDMP holds for less 
restrictive assumptions: t —>■ F{x,t) can be assumed continuous rather than Lipschitz as 
stated above. The same trick also allows to deal with the case where Rtot is time-dependent 
as well. 


( 2 .. 


3. Theoretical result. A general simulation algorithm, called the True Jump Method 
flGT13l . is provided in appendix [A] it is also a general method to prove existence of PDMPs 
flGT13l under the conditions stated in the previous section. A serious limitation to this algo¬ 
rithm is the need to solve the integral equation in step 3. The limitation pertains to the cas^ 
T = 0 if the total transition rate function Rtot is time-dependent. In general, even if the 
flow (</>*) is known analytically, the integral equation needs to be solved numerically using 
Newton method or the bisection algorithm for example. Here, we propose a solution based 
on an ODE solver and a change of time variable. 

Theorem 3.1. Let us make the general assumptions: 

• the flow (()>*)t>o is globally defined, 

• Rtot{x) > 0 for all x € 

• X —t Rtot{x) is locally Lipschitz on 

Then, the solution {X{T ~), T„) of the integral problem in step 3. is given by {y{Sn), t(S'„)) 
where 


y{s) = F{y{s)) / Rtot{y{s)) 

t{s) = 1 / Rtot{y{s)) (3.1) 

t/(0) = X(T„_i), r(0)=T„_i. 


^the continuous function is stationary in this case 
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Also X{t{s)) = y{s) for s e [0, 5'„). 

Proof. The trick for solving jj, ^ Rtot{X{s))ds = 5'„, in Algorithm 1., consists in the 
introduction of the variable r(s) > T„_i defined as 



Rtot{X{t))dt 


= s 


which gives the second equation in ( |3.1| i by differentiating with respect to s. However, we 

ds f 

need to know X{t{s)) to solve the ODE for t(s). Writing y{s) = X{t{s)) and using chain 
rule differentiation gives the remaining equation in (|3.1|l. □ 


The fact that the flow f is globally defined (no ’’explosions”) is a simplifying assumption 
of IIDav84ll . The theorem states that solving the integral equation in step 3. of the Algorithm 
1. amounts to integrate an ODE over a time interval [0, Sn) where Sn is randomly drawn 
from an exponential distribution 8{1) at each recursion of the algorithm. 

Remark 2. It is easy to generalize the theorem to the case where F, Rtot cire time- 
dependent. One finds 


f y{s) = F{y{s),T{s)) / Rtot{y{s), t(s)) 
1 'j'(s) = 1 / Rtot{yis),T{s)) 


3.1. Sampling the continuous flow. The previous method is effective for computing the 

jump instants and the values of the PDMP at these instants. However, if one seeks to sample 
the continuous part in between jumps, there is a difficulty because one only has access to X 
through a change of time variable r. A cheap way around this issue is by adding a Poisson 
process of constant rate A > 0 to our PDMP, thereby increasing the dimension 

of our system: d ^ d 1. This way. Theorem 1. condition Rtot{x) > 0 is always satisfied 
because Rtot{x) > A > 0. This additional jump process will add jump events at rate A, hence 
sampling the continuous dynamics. 

3.2. Simulations. The simulation of the process {X{t)) consists in applying the Algo¬ 
rithm 1. together with the use of Theorem 1. We call this procedure the CHV method. The 
performance of the ODE solver strongly impacts the speed and ’’precision” of the simulation. 
It is expected that the equations eB may be stiff depending on whether a lot (resp. a few) 
reactions are initiated in which case the factor 1 / Rtot might assume vanishing or extremely 
large values. Hence, we find it convenient to use the ODE solvers IIHin831 for their abil¬ 
ity to automatically switch between stiff and non-stiff methods. Note that depending on the 
problem, other ODE methods IIHW96II may be more appropriate than the one suggested. 

3.3. Comparison with the FJM. As stated in the introduction, when the flow is known 
analytically, there is no point in choosing a quasi-exact method e.g. the True Jump Method 
instead of the exact EJM, unless of course, if the jump instants are known analytically. 

When the flow is not known analytically, both methods EJM/CHV are quasi-exact in 
that they bear a numerical error coming from integrating the ODE underlying the continuous 
component. The question we ask is which method is the fastest given a numerical tolerance. 

Eor convenience, we recall the two ODEs that we need to solve between the jumps, on 
an interval of length Sn ~ ^(1). for each method: 


(FJM) : X = F{x)/X or {CHV) : 


y = F{y)/Rtot{y) 
T = l/Rtot{y) 


(3.3) 
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If A « Rtot{x), then (FJM) is smaller hence faster to solve. However, if Rtot{x) varies a 
lot, it may be difficult to bound it effectively using a single constant A as required in the FJM 
(see appendix [B). Hence, the choice of the method boils down to know whether it is faster 
(and more accurate) to solve k time^(FJM) or a single time (CHV). This depends a lot on 
the time stepper and the way it copes with the numerical tolerances to produce the (adaptive) 
time steps. 

In the case of fixed time steps, one needs to see whether the cost of the larger system 
(CHV) compensates the cost of the generation of the k random numbers. Hence, if the FJM 
does not require a lot of fictitious jumps, this method is more appropriate than CHV. 

Finally, if the transition rate Rtot is zero on a open set, then the CHV method breaks 
down and the FJM is more appropriate. On a practical side, the CHV method is simpler to 
program but the FJM assume far less regularity for Rtot- 

4. Numerical examples. We give two examples where our new method is particularly 
effective compared to the standard use of event location for solving the integral equation to 
compute the jumping times. Then we provide an example to compare our method to the FJM. 

4.1. Stiff problem with possibly infinite jumping times. The purpose of this example 
is to show a case where the event location in Mat lab is actually doing a good job and 
compare it to our method. To make things difficult for the event location, we chose stiff 
dynamics between the jumps. Note that the next jumping time can be infinite in our example. 
We look at the following process of switching dynamics where X{t) = {xc{t),Xd{t)) G 

/ ic = a{t)xc .... 

\ a{t) = —100 X {2{xd{t) mod 2) — 1) 


starting from a;c(0) = 1. The (unique) transition rate is Rtot{xc:Xd) = Xc associated with 
the jump Xd ^ Xd + I for which T\{{xc,Xd),dy) = (dy). As Xc increases, the 

transition rate also increases. Hence, the continuous problem is stiff because one has to make 
sure (numerically) that R> 0 despite the abrupt decrease of Xc- Also, the integral equation is 
difficult to handle numerically with high precision. The jumping time, solution of the integral 
equation, reads (see Theorem 3.11; 


def 


a = —100 X {xd{Tn-i) mod 2 — 1) 

1 


Tn = Tn-i + - log (1 + aS'„/xc(T„_i)). 


X(r„) =X(r„_i)+a5„. (4.2) 

From the formula, it appears that the next jumping time can be infinite when l+aiS'„/xc(T„_i) < 


0 when a < 0. The next Figure 4.1 shows an example of realization where the first 15000 
jumping times are computed. To compare the methods, we draw 15000 realizations of a 
random variable with exponential distribution and use this sequence for the 4 examples. 

Note that we had to find a long enough trajectory because most ends quickly as the 
jumping time can be infinite. We see that for similar tolerances, our new method, called 
CHV, performs better as the number of jump events increases. Surprisingly, ode4 5, which is 
not designed for stiff systems performs quite well. We tried all available steppers in Matlab 
and only plot the two best results. 

The following table |4.1| provides the rvinning times (ratios) for the cases in Figure [4~T] 
Despite the additional dynamical variable r of our method, the fastest running times using 
either methods are almost the same. 


'*for k fictitious jumps 
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Fig. 4.1. Top: example of trajectory plotted using analytical formula ^4.2) . Bottom are errors in jumping 
time and amplitude of the process using the event location capabilities of Mat lab (termed event here) and using 
our method (termed CHV). The horizontal axis is the jump index. In both cases, the 2 smallest errors were obtained 
using the steppers odell3 and ode45. The absolute/relative tolerance are atol = le — 12 and rtol = le — 12. 



ode45 

odell3 

event 

1 

1.851 

CHV 

1.009 

1.0967 


Table 4.1 

Computation time ratios for the examples in Figure^j\ 


4.2. Stiff problem with finite jumping times. We now provide a numerical example to 
show that the scheme in IIRiel3l can lead to wrong dynamics. We chose the following process 
of switching dynamics where X{t) = {xc{t),Xd{t)) G 


ic = 10 • X, 

Xc = —‘ixl 


if 

otherwise 


Xd even 


(4.3) 
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where the (unique) transition rate is Rtot{xc,Xd) = Xc associated with the jump Xd —>■ 
ccrf + l forwhichn((a;c,a:£;),(ij/) = < 5 ( 2 ;^ 3 ,^+ 1 ) (<^2/)- The jumping time solution of the integral 
equation reads (see Theorem |3.1|i: 


'X{Tr,) 

ATn), 


10 • <S'„ + Xc{Tn-l) 

10S„ 


ra log (1 + 




.3a;(T„_i) 


(e3S. _ 


if Xd even 


otherwise. 


(4.4) 


We see that the function t is well-defined meaning that the jumping times are always finite. 
The results are plotted in Figure [4~2| To compare the methods, we draw 20 realizations of a 


Process trajectory 



Jumping time error (event) 



Amplitude error (event) 



Jumping time error (CHV) 



Amplitude error (CHV) 



Fig. 4.2. Top: example of trajectory plotted using analytical formula \4.4\ . The green square are the jumps 
computed with the event function, the red crosses with our method. Bottom are errors in jumping time and amplitude 
of the process using the event location capabilities of Matlab (termed event here) and using our method (termed 
CHV). The horizontal axis is the jump index. In both cases, the 2 smallest errors were obtained using the steppers 
odell3 and ode45. The absolute/relative tolerance are atol = le — 12 and rtol = le — 12. 


random variable of exponential distribution and use this sequence for the 4 examples. After 
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a few jumps (on the order of ten), the method using the Event function in Matlab fails to 
compute accurately a jumping time when Xc is close to zero. This difficulty in computing 
the jumping time is also reflected in the cpu time shown in Table |4.2] Hence, computing 100 
jumps with the Event method leads to completely wrong dynamics whereas our new method 
gives very satisfactory results when compared to the analytical results. 



ode45 

odell3 

event 

4.941 

3.569 

CHV 

1.940 

1 


Table 4.2 

Computation time ratios for the examples in figure [42] 


4.3. Comparison with the rejection method. Finally, we provide a numerical example 
to compare our method to the FJM. We chose the following process of switching dynamics: 

{ ±c = 3 ■ Xc if Xd even 

Xc = —4:Xc otherwise (4.5) 

:cc(0) = 0.05, a:d(0) = 0 

where the (unique) transition rate is i?tot(a;c,a:d) = o_)_gi,,,+5.o +0.1 with n((a;e, Xd), dy) = 
5{xa,xi+i){dy)- When the continuous component Xc is increasing, the total rate Rtot is 
bounded by 1. However, when it is decreasing, the total rate is bounded by Rtot{xc{Tn) ^ XdiTn)) 
for t > Tn where T„ is the last jumping time. Hence, we expect the FJM to perform less fic¬ 
titious jump when Xd is odd. The fact that Xc grows exponentially disfavors the FJM as we 
shall see. 

To compare the two methods CHV and FJM, we use the adaptive ODE solver CVODE 
from SUNDIAFS llHBG~*~05ll for its ability to automatically switch between stiff and non¬ 
stiff methods. In each of the 4 cases, we compute 1 jump instant (resp. 3 jump instants) with 
each method and we perform 10® realizations. We show the results in Figure |43] where each 
plot is an histogram of the computing time in seconds of the jump instants. A first interesting 
result concerns the case FJM - 1 jump where it seems that we can see, in the distribution, 
the fictitious jumps used to compute the first jump. Indeed, in our implementation, it takes 
roughly 1.5tos to solve the ODE with both methods. The top two plots show that the CHV 
method is faster for the reasons stated at the beginning of the section. The other 2 plots show 
that this result is not affected by the computation of the jumping time in the case Xd odd, 
where the FJM should be slightly faster. Hence, the fastest method for simulating the process 
( |4.5| l would be to use our method when Xd is even and the rejection method otherwise. 

5. Generalizations. Several possible generalizations of our model suggest themselves. 
Basically, these entail a modification of the flow ( |2.2] l or of the jump process. 

If the flow can be written as a Cauchy problem, then our method still applies albeit 
possibly in a different state space. For example, we can consider the case of a flow generated 
by a delay differential equation IIHal93l or if jump process depends on the history of the 
process {X{t)). These cases are handled very similarly to the case presented here, except 
that \2.2\ , ( |3.1| l are now delay differential equations. One could also look at the case of a flow 
generated by partial differential differential equations IIBRl II . 

The next possible generalization was introduced in IIBVTH05I where a time delay was 
added between the initiation and the completion of some, or all, of the reactions. These delays 
were motivated by oscillations observed in gene regulation networks BDFDOlll . A modified 
True Jump Method was given IIBVTHOSI where a table is added to record when the delayed 
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Fig. 4.3. Histograms of the computing time in seconds for each method CHV or FJM when computing 1 jump 
(top) and 3 jumps (bottom). The absolute/relative tolerance are atol = le — 10 and rtol = le — 10. See text for a 
description of the stochastic process. 


reactions are scheduled. The algorithm was later improved by OCai071lAnd07l but the method 
is the same. Hence, the computation of the next reaction time follows from the same integral 
equation that we solved here and our results equally apply. 

6. Discussion. We have identified a new way to simulate a large class of Markov pro¬ 
cesses, the PDMPs, with a method which is versatile enough to be applied to various con¬ 
tinuous part dynamics (ODE / PDE / delay differential equations) and any continuous rate 
function Rtot- Indeed the precision of the method is directly linked to the one of the ODE 
solver, which we assume is arbitrary given the state of the art concerning global error estima¬ 
tion OSkeShinOsOTl . 

As a particular case, our method provides a solution to cope with the notoriously difficult 
problem IIAnd07l of the simulation of chemical reactions with time dependent propensity 
functions. Nevertheless, this problem is only formal as only the Eictitious Jump Method is 
exact in this case. 

Our result is an important achievement for the following reasons. Eirst, there is no need 
to modify existing ODE solvers to implement our method. Second, it provides a way to use 
an adaptive discretization to approximate both the ODE flow and the integral equation. As 
such, it makes the method llACT~*~05l IRiel3l for the simulation of PDMPs not relevant. It 
also allows to focus on the choice of the ODE solver without ad-hoc modification of existing 
tools, this is especially useful when studying stiff systems such as the slow-fast ones, largely 
present in the modeling of biology. Our results also makes it easy to consider flows driven by 
PDE or delay differential equations for example, without sacrificing for the numerical preci- 
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sion. Further, our results bridge the gap between the numerical methods for the simulation of 
PDMPs and the ones of ODE, making the hrst a subproblem of the second class. In particular, 
numerical error analysis IIHW961ISke861l^s01IICP04l of ODE directly maps onto PDMPs. 

The only relevant question which remains is whether to chose our method or the EJM 
when the continuous flow is not known analytically and when the total rate function Rtot is 
continuous. This is strongly dependent on the total rate function variations and not so much 
on its absolute values. If the probability of having at least 2 fictitious jumps is ’high’, then 
our new method might be the most efficient while being simpler to program. Einally, note 
that mixing the two methods can be effective by choosing the fastest one depending on the 
value of the discrete component as we suggested in our third example. 
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Appendix A. True Jump Method. We recall the true jump method 0GT13I for the exact 
simulation of {X (<))tg]i{+ on a time interval [0, T], This is the SSA Doob-Gillespie algorithm 
IIGil77l when F = 0, for which cj)* = Id. 

Algorithm 1. (True Jump Method) 

1. Draw the initial value Xq 

2. Draw S'„ according to an exponential distribution £{1) 

3. Setr„ = inf{M>T„_i : andX(f) = (^‘(A(T„_i)) 

for T„_i <t <Tn. 

4. Draw X{Tn) according to the law (A(T„_i)), dy). 

5. If T„ < T, return to step 2. 


Appendix B. Fictitious Jump Method. We recall the ’’rejection” method IIGT13I for 
the exact simulation of {X{t))t^^+ on a time interval [0, T], We assume that sup Rtot{x) < 

a;GR‘^ 

A < oo. 

Algorithm 2. (Eictitious Jump Method) 

1. Draw the initial value Xq 

2. Draw S'„ according to an exponential distribution £(A) 

3. Set X{t) = <^‘(A(r„_i)) for r„_i < t < r„_i + Sn. 

4. • either, with probability 1 - , set X(T^) = X(T^-) 

• or else draw X{Tn) according to the law n(A(r„—), dy) 

5. If T„ < T, return to step 2. 
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